# include "scene.h"
# include <stdio.h>
# include <malloc.h>
# include <string.h>
# include <math.h>
# include <errno.h>

# ifdef _OPENMP_
# include <omp.h>
# endif

extern short conf_grid(struct grid_config *gc, struct grid *g);
extern short grid_gen(struct grid_config *gc, struct grid *g);

# ifdef SOLVE_ENERGY
extern short heater_param(struct grid *g);
# endif

static short print_grid_conf(const struct grid_config *gc, const struct grid *g);

extern void mem_alloc(FILE *lptr, struct arr_struct *arr);
extern void mem_distroy(FILE *lptr, struct arr_struct *arr);

extern short solve(struct grid *g);

int main(int argc, char *argv[])
  {
   struct grid g;
   struct grid_config gc;
   struct arr_struct tnext;

   printf("Program starts: in 'main()'...\n");

   g.lptr = stderr;

   if( (g.lptr = (LOG_FILE ? fopen(LOG_FILE, "a"):stderr)) == (FILE *)0) g.lptr = stderr;
   fprintf(g.lptr, "\n\n|----------[%d]. PROGRAM STARTS FROM 't = %E Sec'---------|\n\n", RUN_CNT, t_INIT);

#  ifdef TRACE_FILE
   if( (g.trptr = (TRACE_FILE ? fopen(TRACE_FILE, "w"):stderr)) == (FILE *)0 )
       {
        fprintf(g.lptr, "ERROR:: In 'main()': Unable to open '%s' for writing: %s\n", TRACE_FILE, strerror(errno));
        return -1;
       }
#  endif

   g.nCPU = 1;

#  ifdef _OPENMP_
   g.nCPU = omp_get_num_procs();
   omp_set_num_threads(g.nCPU);
   fprintf(g.lptr, "Parallel computation using OpenMP.\n");
#  endif
   fprintf(g.lptr, "Number of CPUs available for the OpenMP = %d\n", g.nCPU);

/* --- Grid generation--- */
   if(conf_grid(&gc, &g) < 0) return -1;

   if(print_grid_conf((const struct grid_config *)&gc, (const struct grid *)&g) < 0) return -1;

/* Allocating arrays to store the grid-POINTS. */
   g.z = (double *)malloc((g.nz + 2*NzG + 1)*sizeof(double));
   g.r = (double *)malloc((g.nr + 2*NrG + 1)*sizeof(double));
   g.th = (double *)malloc((g.nth + 2*NthG + 1)*sizeof(double));

/* Files to store grid-POINTS. */
   gc.z_file = Z_FILE;
   gc.r_file = R_FILE;
   gc.th_file = TH_FILE;

   if(grid_gen(&gc, &g) < 0) return -1;

/* --- Grid generation done --- */

#  ifdef SOLVE_ENERGY
   fprintf(g.lptr, "NOTE: The energy equation will also be solved. Heater will be ON after %E seconds of the flow-start.\n", t_HEATER_ON);

/* Finding the heater's radius on the top surface. */
   if(heater_param(&g) < 0) return -1;
#  endif


#  ifdef CHECK_CONFIG
   printf("\nATTENTION:\n\t1. t_INIT = %E Sec\n", t_INIT);

// #    ifdef DELTAt
     printf("\t2. DELTAt = %E Sec\n", DELTAt);
// #    else
//      printf("\t2. CLF_INIT = %E Sec\n", CFL_INIT);
// #    endif

#    ifdef RECOV_FILE
     printf("\t3. Make sure that the recovery file, '%s', corresponds to the time-instant '%E Sec'.\n", RECOV_FILE, t_INIT);
     if(!strcmp(RECOV_FILE, W_RECOV_FILE)) printf("\t4. RECOV_FILE = W_RECOV_FILE = '%s'. Recommended: RECOV_FILE != W_RECOV_FILE.\n", RECOV_FILE);
#    endif

   printf("\nAre all these parameters reasonable? Continue (y/n) ? "); fflush(stdout);
   if(getchar() != 'y') return 0;
   else while(getchar() != '\n');
#  endif

/* Allocating memory. */

/* Allocating memory to store 3xp and 4x(u, v, w), and 1xT, while solving the Momentum and Energy equations. */
   g.arr.n1 = g.nr + 2*NrG;
   g.arr.n2 = g.nth + 2*NthG;
   g.arr.n3 = g.nz + 2*NzG;
   g.arr.nq = 15;
# ifdef SOLVE_ENERGY
   g.arr.nq = 16;
# endif
   mem_alloc(g.lptr, &g.arr);

/* Allocating memory to store A at the cell-faces*/
   g.memFace_r.n1 = g.nr + 1;
   g.memFace_r.n2 = g.nth;
   g.memFace_r.n3 = g.nz;
   g.memFace_r.nq = nRK;
   mem_alloc(g.lptr, &g.memFace_r);

/* Allocating memory to store B at the cell-faces*/
   g.memFace_th.n1 = g.nr;
   g.memFace_th.n2 = g.nth+1;
   g.memFace_th.n3 = g.nz;
   g.memFace_th.nq = nRK;
   mem_alloc(g.lptr, &g.memFace_th);

/* Allocating memory to store C at the cell-faces*/
   g.memFace_z.n1 = g.nr;
   g.memFace_z.n2 = g.nth;
   g.memFace_z.n3 = g.nz+1;
   g.memFace_z.nq = nRK;
   mem_alloc(g.lptr, &g.memFace_z);

/* Allocating memory to store u,v,w,T at the outlet for t_BACKUP number of time-steps.
   It is required for application of the Orlansky BC. */
   g.memOut.n1 = g.nr_gap;
   g.memOut.n2 = g.nth;
   g.memOut.n3 = zCells_BACKUP;
   g.memOut.nq = (DEPVAR_N-1)*t_BACKUP;
   mem_alloc(g.lptr, &g.memOut);

/* Allocating memory to store dependent variables at outlet-cell. It will be used in application of the Orlansky BC. */
   tnext.n1 = DEPVAR_N-1; /* Corresponding to u, v, w, T */
   tnext.n2 = g.nr_gap;
   tnext.n3 = g.nth;
   tnext.nq = 1;
   mem_alloc(g.lptr, &tnext);
   g.outBC.tnext_ug = tnext.q[0][0];
   g.outBC.tnext_vg = tnext.q[0][1];
   g.outBC.tnext_wg = tnext.q[0][2];
#  ifdef SOLVE_ENERGY
   g.outBC.tnext_Tg = tnext.q[0][3];
#  endif

/* Memory allocation done. */

   solve(&g);

/* Freeing memory. */
   mem_distroy(g.lptr, &g.arr);
   mem_distroy(g.lptr, &g.memFace_r);
   mem_distroy(g.lptr, &g.memFace_th);
   mem_distroy(g.lptr, &g.memFace_z);
   mem_distroy(g.lptr, &g.memOut);
   mem_distroy(g.lptr, &tnext);
//    mem_distroy(g.lptr, &coeff);


#  ifdef TRACE_FILE
   if(g.trptr != stderr) fclose(g.trptr);
#  endif

   fprintf(g.lptr, "...'main()' ends!\n");
   if(g.lptr != stderr) fclose(g.lptr);

   printf("...program ends.\n");

   return 0;
  }


static short print_grid_conf(const struct grid_config *gc, const struct grid *g)
  {
   FILE *lptr;

#  ifdef TRACE_FILE
   fprintf(g->trptr, "Entering into 'print_grid_conf()'...\n");
#  endif

   lptr = g->lptr;


#ifdef CHECK_CONFIG
top:
#endif

   fprintf(lptr, "\n\n['r-grid' parameters]\n\n");

#  ifdef r_UNIF_GRID
   fprintf(lptr, "dr: %E ---> %E\n", dR_CYL, g->dr);
#  else
   fprintf(lptr, "dr_center = %E\n", dr_CENTER);
   fprintf(lptr, "r_cd_center ---> %E\n", gc->r_cd_center);
   fprintf(lptr, "**nr_center ---> %d\n", gc->nr_center);
   fprintf(lptr, "*r_mid1: %E ---> %E\n", r_MID1, gc->r_mid1);
   fprintf(lptr, "dr_mid1: %E ---> %E\n", dr_MID1, gc->dr_mid1);
   fprintf(lptr, "dr_mid2 = %E\n", dr_MID2);
   fprintf(lptr, "*r_mid2: %E ---> %E\n", r_MID2, gc->r_mid2);
   fprintf(lptr, "r_cd_mid12 ---> %E\n", gc->r_cd_mid12);
   fprintf(lptr, "**nr_mid12 ---> %d\n", gc->nr_mid12);
   fprintf(lptr, "dr_base: %E ---> %E\n", dr_BASE, gc->dr_base);
   fprintf(lptr, "r_cd_base ---> %E\n", gc->r_cd_base);
   fprintf(lptr, "**nr_base ---> %d\n", gc->nr_base);
   fprintf(lptr, "*r_base = %E\n", r_BASE);
   fprintf(lptr, "dr_wall: %E ---> %E\n", dr_WALL, gc->dr_wall);
   fprintf(lptr, "r_cd_wall ---> %E\n", gc->r_cd_gap);
   fprintf(lptr, "**nr_gap ---> %d\n", g->nr_gap);

#  endif
   fprintf(lptr, "*r_cylinder = %E\n", r_CYL);

   fprintf(lptr, "\n**nr_nozzle ---> %d\n", g->nr_nozzle);
   fprintf(lptr, "**nr_solid ---> %d\n", g->nr_solid);
   fprintf(lptr, "**nr_gap ---> %d\n", g->nr_gap);
   fprintf(lptr, "\n**nr ---> %d\n", g->nr);


   fprintf(lptr, "\n\n['th-grid' parameters]\n\n");

#  ifdef th_UNIF_GRID
   fprintf(lptr, "d_theta:\n\t%E ---> %E (in degrees),\n\t%E ---> %E (in radians)\n",\
                  dTHETA, 90*g->dth/acos(0), dTHETA*acos(0)/90, g->dth);

   fprintf(lptr, "**nth = %d\n", g->nth);
#  else
   fprintf(lptr, "ERROR:: in 'print_grid_conf()': incomplete definition for non-uniform grids in thera direction. Exiting the program.\n");
   exit(1);
#  endif


   fprintf(lptr, "\n\n['z-grid' parameters]\n\n");

   fprintf(lptr, "h_cylinder = %E\n", h_CYL);

#  ifdef z_UNIF_GRID
   fprintf(lptr, "dz: %E ---> %E\n", dZ, g->dz);
#  else
   fprintf(lptr, "dz_top = %E\n", dz_TOP);
   fprintf(lptr, "**nz_top ---> %d\n", gc->nz_top);
   fprintf(lptr, "z_cd_top ---> %E\n", gc->z_cd_top);
   fprintf(lptr, "dz_mid2 = %E\n", dz_MID2);
   fprintf(lptr, "*z_mid2: %E ---> %E\n", z_MID2, gc->z_mid2);
   fprintf(lptr, "**nz_mid12 ---> %d\n", gc->nz_mid12);
   fprintf(lptr, "z_cd_mid12 ---> %E\n", gc->z_cd_mid12);
   fprintf(lptr, "dz_mid1 = %E\n", dz_MID1);
   fprintf(lptr, "*z_mid1: %E ---> %E\n", z_MID1, gc->z_mid1);
   fprintf(lptr, "**nz_mid1 ---> %d\n", gc->nz_mid1);
   fprintf(lptr, "z_cd_bot ---> %E\n", gc->z_cd_bot);
   fprintf(lptr, "**nz_unif_bot = %d\n", nz_UNIF_BOT);
   fprintf(lptr, "dz_unif_bot: %E ---> %E\n", dz_UNIF_BOT, gc->dz_unif_bot);
#  endif
   fprintf(lptr, "\n**nz ---> %d\n\n", g->nz);

   fflush(lptr);


#  ifdef CHECK_CONFIG
   if(lptr != stdout)
       {
        lptr = stdout;
        goto top;
       }
   else{
        fprintf(lptr, "\tAre the parameters reasonable (y/n) ? ");  fflush(lptr);
        if(getchar() != 'y')
           {
            fprintf(lptr, "\nWARNING:: in 'print_grid_conf()': Re-define the grid-configuration parameters in \"scene.h\", re-compine, and then run the program again.\n\n");

            while(getchar() != '\n');
            return -1;
           }
        else while(getchar() != '\n');

        lptr = g->lptr;
       }
#  endif

#  ifdef TRACE_FILE
   fprintf(g->trptr, "...'print_grid_conf()' ends.\n");
#  endif

   return 1;
  }
